# Import PyDicom before using PyDicom functions
# DICOM (Digital Imaging and Communications in Medicine) is the standard protocol for
# the management and transmission of medical images and related data, and is used by
# many healthcare facilities
# PyDicom is a pure Python package for working with DICOM files
import pydicom
# The `pathlib` module is similar to the `os.path` module, but `pathlib` provides a more
# advanced and convenient interface than `os.path`
# It is possible to use `pathlib` to represent file paths as specialized `Path` objects
# instead of plain strings
from pathlib import Path
# SimpleITK is a simplified programming interface to the algorithms and data structures
# of the Insight Toolkit (ITK) for segmentation, registration and advanced image analysis
# SimpleITK supports bindings for multiple programming languages including C++, Python, R,
# Java, C#, Lua, Ruby and TCL
# Combining SimpleITK’s Python bindings with the Jupyter notebook web application creates
# an environment which facilitates collaborative development of biomedical image analysis
# workflows
# In this project, SimpleITK provides the ability to automatically detect and read all DICOM
# files so that the user does not have to manage file reading or slice sorting
import SimpleITK as sitk
import numpy as np
import torch
from torchvision.utils import make_grid
import matplotlib as mpl
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
from datetime import datetime
from functools import wraps
import itertools
import math
import random
import reprlib
import sys
%matplotlib inline
XINHUI = "#7a7374"
XUEBAI = "#fffef9"
YINBAI = "#f1f0ed"
YINHUI = "#918072"
figure_size = (16, 9)
custom_params = {
"axes.axisbelow": True,
"axes.edgecolor": YINBAI,
"axes.facecolor": XUEBAI,
"axes.grid": True,
"axes.labelcolor": XINHUI,
"axes.spines.right": False,
"axes.spines.top": False,
"axes.titlecolor": XINHUI,
"figure.edgecolor": YINBAI,
"figure.facecolor": XUEBAI,
"grid.alpha": 0.8,
"grid.color": YINBAI,
"grid.linestyle": "--",
"grid.linewidth": 1.2,
"legend.edgecolor": YINHUI,
"patch.edgecolor": XUEBAI,
"patch.force_edgecolor": True,
"text.color": XINHUI,
"xtick.color": YINHUI,
"ytick.color": YINHUI,
}
mpl.rcParams.update(custom_params)
reprlib_rules = reprlib.Repr()
reprlib_rules.maxother = 250
sys.path.append("../")
from Modules import *
# The dataset used in this practice project is a very small subset of CT images extracted from
# the Cancer Imaging Archive (TCIA), which contains intermediate slices of all CT images
# in which valid age, mode, and contrast markers can be found
# The dataset is provided by the Kaggel Datasets called CT Medical Imaging
# (https://www.kaggle.com/datasets/kmader/siim-medical-images),
# license type is Attribution 3.0 Unported (CC BY 3.0)
# https://creativecommons.org/licenses/by/3.0
# The link to the TCIA archive of the full dataset is
# https://wiki.cancerimagingarchive.net/display/Public/TCGA-LUAD
dir_path = "../Datasets/Kaggle - CT Medical Images/dicom_dir/"
sample_dcm = "ID_0000_AGE_0060_CONTRAST_1_CT.dcm"
# The main function of PyDicom to read and parse DICOM files is `read_file`
dicom_file = pydicom.read_file(dir_path + sample_dcm)
tabulation = Form_Generator()
tabulation.heading_printer("Initial understanding of DICOM file")
statements = [
"""
dir_path = "../Datasets/Kaggle - CT Medical Images/dicom_dir/"
sample_dcm = "ID_0000_AGE_0060_CONTRAST_1_CT.dcm"
dicom_file = pydicom.read_file(dir_path + sample_dcm)
"""
]
tabulation.statement_generator(statements)
variables = ["dicom_file"]
values = [str(reprlib_rules.repr(dicom_file))]
tabulation.variable_generator(variables, values)
expressions = [
"dicom_file[0x0028, 0x0010]",
"dicom_file[0x0028, 0x0011]",
"dicom_file[0x0018, 0x0015]",
"dicom_file.Rows",
"dicom_file.Columns",
"dicom_file.BodyPartExamined",
"dicom_file.keys()",
"dicom_file.values()",
"dicom_file.dir()",
'dicom_file.dir("Image")',
]
results = [
str(dicom_file[0x0028, 0x0010]),
str(dicom_file[0x0028, 0x0011]),
str(dicom_file[0x0018, 0x0015]),
str(dicom_file.Rows),
str(dicom_file.Columns),
str(dicom_file.BodyPartExamined),
str(reprlib_rules.repr(dicom_file.keys())),
str(reprlib_rules.repr(dicom_file.values())),
str(reprlib_rules.repr(dicom_file.dir())),
str(reprlib_rules.repr(dicom_file.dir("Image"))),
]
tabulation.expression_generator(expressions, results, 1)
Initial understanding of DICOM file +-------------------------------------------------------+ | Statement | +-------------------------------------------------------+ | dir_path = "../Datasets/Kaggle - CT Medical | | Images/dicom_dir/" | | sample_dcm = "ID_0000_AGE_0060_CONTRAST_1_CT.dcm" | | dicom_file = pydicom.read_file(dir_path + sample_dcm) | +-------------------------------------------------------+ +------------+------------------------------------------------+ | Variable | Value | +------------+------------------------------------------------+ | dicom_file | Dataset.file_meta | | | ------------------------------- | | | (0002, 0000) File Meta Information Group | | | Length UL: 194 | | | (0002, 0001) Fil...Group Length | | | UL: 524296 | | | (7fe0, 0010) Pixel Data | | | OW: Array of 524288 elements | +------------+------------------------------------------------+ +-----------------------------+-------------------------------+ | Expression | Result | +-----------------------------+-------------------------------+ | dicom_file[0x0028, 0x0010] | (0028, 0010) Rows | | | US: 512 | | dicom_file[0x0028, 0x0011] | (0028, 0011) Columns | | | US: 512 | | dicom_file[0x0018, 0x0015] | (0018, 0015) Body Part | | | Examined | | | CS: 'CHEST' | | dicom_file.Rows | 512 | | dicom_file.Columns | 512 | | dicom_file.BodyPartExamined | CHEST | | dicom_file.keys() | dict_keys([(0008, 0000), | | | (0008, 0005), (0008, 0008), | | | (0008, 0016), (0008, 0018), | | | (0008, 0020), (0008, 0021), | | | (0008, 0022), ...095, 0010), | | | (0097, 0000), (0097, 0010), | | | (0099, 0000), (0099, 0010), | | | (7003, 0000), (7003, 0010), | | | (7fe0, 0000), (7fe0, 0010)]) | | dicom_file.values() | dict_values([(0008, 0000) | | | Group Length | | | UL: 430, (0008, 0005) | | | Specific Character Set | | | CS:...up Length | | | UL: 524296, | | | (7fe0, 0010) Pixel Data | | | OW: | | | Array of 524288 elements]) | | dicom_file.dir() | ['AccessionNumber', | | | 'AcquisitionDate', | | | 'AcquisitionNumber', | | | 'AcquisitionTime', | | | 'BitsAllocated', | | | 'BitsStored', ...] | | dicom_file.dir("Image") | ['ImageDimensions', | | | 'ImageFormat', | | | 'ImageGeometryType', | | | 'ImageLocation', | | | 'ImageOrientation', | | | 'ImageOrientationPatient', | | | ...] | +-----------------------------+-------------------------------+
# DICOM data has an attribute called `pixel_array` that provides more useful pixel data
# for uncompressed images that can be passed to the graphics library for viewing
# To use this attribute, the system must have the NumPy numeric package installed,
# since `pixel_array` returns a NumPy array
ct = dicom_file.pixel_array
tabulation = Form_Generator()
tabulation.heading_printer("Getting pixel data from DICOM file")
statements = ["ct = dicom_file.pixel_array"]
tabulation.statement_generator(statements)
variables = ["ct"]
values = [str(reprlib_rules.repr(ct))]
tabulation.variable_generator(variables, values)
expressions = ["ct.shape"]
results = [str(ct.shape)]
tabulation.expression_generator(expressions, results)
Getting pixel data from DICOM file +-----------------------------+ | Statement | +-----------------------------+ | ct = dicom_file.pixel_array | +-----------------------------+ +----------+------------------------------------------------+ | Variable | Value | +----------+------------------------------------------------+ | ct | array([[0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | ..., | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0]], dtype=uint16) | +----------+------------------------------------------------+ +------------+------------+ | Expression | Result | +------------+------------+ | ct.shape | (512, 512) | +------------+------------+
def image_display(image, ax, title, cmap):
ax.imshow(image, cmap)
ax.grid(False)
ax.set_title(title, loc="center", pad=10)
x_ticks = list(range(0, image.shape[1], 50))
y_ticks = list(range(0, image.shape[0], 50))
ax.set(xticks=x_ticks, xticklabels=x_ticks,
yticks=y_ticks, yticklabels=y_ticks)
ax.set_xlim(left=0)
ax.set_ylim(top=0)
ax.minorticks_on()
return ax
# Path classes are divided between pure paths, which provide purely computational
# operations without I/O, and concrete paths, which inherit from pure paths but also
# provide I/O operations
path_object = Path(dir_path)
# `PurePath.name` returns a string representing the final path component, excluding
# the drive and root directory (if any)
# When a path points to a directory, `Path.iterdir` generates a path object of the directory's
# contents
# `Path.is_file` returns True if the path points to a normal file or to a symbolic link
# to a normal file, False if it points to another type of file
random_dicom_path = random.choice(
[
file.name
for file in path_object.iterdir()
if file.is_file() & (pydicom.read_file(file).BodyPartExamined != "CHEST")
]
)
random_dicom_file = pydicom.read_file(dir_path + random_dicom_path)
plt.rcParams["figure.figsize"] = (
figure_size[0] / 3 * 2, figure_size[1] / 4 * 5)
fig, axs = plt.subplots(nrows=2, ncols=2)
image_display(
ct,
axs[0, 0],
f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'gray'",
cmap="gray",
)
image_display(
ct,
axs[0, 1],
f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'bone'",
cmap="bone",
)
image_display(
random_dicom_file.pixel_array,
axs[1, 0],
f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'gray'",
cmap="gray",
)
image_display(
random_dicom_file.pixel_array,
axs[1, 1],
f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'bone'",
cmap="bone",
)
fig.suptitle(
"Visual Comparison of CT Medical Image Display Using Different Colormaps",
fontsize="x-large",
x=0.5,
y=0,
)
plt.tight_layout()
plt.show()
def annotation_color(cmap):
cmap_color = mpl.colormaps[cmap]
return cmap_color(0.85)
def first_max(list):
max_element = list[0]
for element in list[1:]:
if len(element) > len(max_element):
max_element = element
return max_element
def text_aligner(text):
line_list = text.split("\n")
max_length = len(first_max(line_list))
line_list = [
(":" + " " * (max_length - len(line))).join(line.split(":"))
if len(line) < max_length
else line
for line in line_list
]
text = "\n".join(line_list)
return text
def plane_information_annotator(func):
@wraps(func)
def wrapper(ax, cmap, dicom_file, fontsize_adjustment):
plane = func(ax, cmap, dicom_file, fontsize_adjustment)
text = f"Image: {dicom_file.Columns} x {dicom_file.Rows}"
text += f"\nImage Plane: {plane}"
text += f"\nPatient Position: {dicom_file.PatientPosition}"
text += f"\nSlice Thickness: {dicom_file.SliceThickness:.1f} mm"
text += f"\nSlice Location: {dicom_file.SliceLocation:.1f} mm"
ax.text(
0.575,
0.825,
text_aligner(text),
transform=ax.transAxes,
fontsize=9 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
return wrapper
def orientation_annotator(func):
@wraps(func)
def wrapper(ax, cmap, dicom_file, fontsize_adjustment):
plane, top, right, bottom, left = func(dicom_file)
ax.text(
0.4875,
0.95,
top,
transform=ax.transAxes,
fontweight="heavy",
fontsize=11 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
ax.text(
0.4875,
0.025,
bottom,
transform=ax.transAxes,
fontweight="heavy",
fontsize=11 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
ax.text(
0.95,
0.4875,
right,
transform=ax.transAxes,
fontweight="heavy",
fontsize=11 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
ax.text(
0.025,
0.4875,
left,
transform=ax.transAxes,
fontweight="heavy",
fontsize=11 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
return plane
return wrapper
def image_orientation(func):
@wraps(func)
def wrapper(*args, **kwargs):
Left, Right, Anterior, Posterior, Head, Feet = "L", "R", "A", "P", "H", "F"
axial_plane = [Anterior, Right, Posterior, Left]
coronal_plane = [Right, Head, Left, Feet]
sagittal_plane = [Anterior, Head, Posterior, Feet]
plane = func(*args, **kwargs)
# Patient Position (0018,5100) specifies the position of the patient relative to
# the space of the imaging device and is used for annotation purposes only,
# but does not provide a precise mathematical relationship between the patient and
# the imaging device
patient_position = dicom_file.PatientPosition
if plane == "Axial":
term_list = ["S", "DR", "P", "DL"]
top = 4 - term_list.index(patient_position[2:])
if patient_position.startswith("HF"):
right = top + 1
left = right + 2
else:
left = top + 1
right = left + 2
bottom = top + 2
top, right, bottom, left = [
axial_plane[index] if index < 4 else axial_plane[index - 4]
for index in [top, right, bottom, left]
]
elif plane == "Coronal":
term_list = ["DL", "DR"]
top = term_list.index(patient_position[2:]) * 2
if patient_position.startswith("AF"):
right = top + 1
left = right + 2
else:
left = top + 1
right = left + 2
bottom = top + 2
top, right, bottom, left = [
coronal_plane[index] if index < 4 else coronal_plane[index - 4]
for index in [top, right, bottom, left]
]
elif plane == "Sagittal":
term_list = ["LF", "RF"]
if patient_position.endswith("S"):
top = 0
right = term_list.index(patient_position[:2]) * 2 + 1
left = right + 2
else:
top = 2
left = term_list.index(patient_position[:2]) * 2 + 1
right = left + 2
bottom = top + 2
top, right, bottom, left = [
sagittal_plane[index] if index < 4 else sagittal_plane[index - 4]
for index in [top, right, bottom, left]
]
else:
top, right, bottom, left = np.repeat("", 4)
return plane, top, right, bottom, left
return wrapper
@plane_information_annotator
@orientation_annotator
@image_orientation
def get_plane_information(dicom_file):
# The DICOM standard defines a patient-oriented reference coordinate system (RCS) that
# enables the user to measure the position and orientation of the image relative to
# the patient
# Image Orientation (Patient) (0020,0037) specifies the cosine of the orientation of
# the first row and column relative to the patient, which should be provided as a pair,
# where the row values for the x, y, and z axes are followed by the column values for
# the x, y, and z axes
# The orientation of the axes is determined solely by the orientation of the patient, i.e.,
# the RCS and Image Orientation (Patient) (0020,0037) values specify the orientation of
# the image frame rows and columns
IOP = dicom_file.ImageOrientationPatient
row_xyz = IOP[:3]
column_xyz = IOP[3:]
# `numpy.cross` returns the cross product of two vector arrays
# The geometric interpretation of the cross product is that two vectors determine a plane,
# and the cross product points in a direction different from both vectors
# The Patient-Based Coordinate System is a right-handed system, i.e., the vector
# cross product of a unit vector along the positive x-axis and a unit vector along
# the positive y-axis is equal to a unit vector along the positive z-axis
plane = [abs(round(x)) for x in np.cross(row_xyz, column_xyz)]
# The Image Type (0008,0008) identifies important image identification features
# and it returns the Pixel Data Characteristics, Patient Examination Characteristics,
# Modality Specific Characteristics, and Implementation Specific Identifiers
# While the third value (Modality Specific Characteristics) should identify any
# specialization specific to the Image Information Object Definition (IOD), it
# and the values that follow it are not mandatory, unlike the first two values,
# and therefore some DICOM data does not have this value
if plane[0] == 1 and plane[1] == 0 and plane[2] == 0:
return "Sagittal"
elif plane[0] == 0 and plane[1] == 1 and plane[2] == 0:
return "Coronal"
elif plane[0] == 0 and plane[1] == 0 and plane[2] == 1:
return "Axial"
else:
"Unknown"
def get_patient_information(ax, cmap, dicom_file, fontsize_adjustment, patien_age):
text = f"Patient ID: {dicom_file.PatientID}"
text += f"\nPatient's Sex: {dicom_file.PatientSex}"
try:
patien_age = dicom_file.PatientAge.strip(str(0))[:-1]
except AttributeError:
patien_age = patien_age
text += f"\nPatient's Age: {patien_age} y"
text += f"\nModality: {dicom_file.Modality}"
text += f"\nBody Part Examined: {dicom_file.BodyPartExamined}"
ax.text(
0.04,
0.825,
text_aligner(text),
transform=ax.transAxes,
fontsize=9 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
def get_date_information(ax, cmap, dicom_file, fontsize_adjustment):
date = datetime.strptime(dicom_file.StudyDate,
"%Y%m%d").strftime("%Y.%m.%d")
text = f"Study Date: {date}"
ax.text(
0.04,
0.05,
text,
transform=ax.transAxes,
fontsize=7 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
def get_manufacturer_information(ax, cmap, dicom_file, fontsize_adjustment):
text = f"Manufacturer: {dicom_file.Manufacturer}"
ax.text(
0.525,
0.05,
text.rjust(36),
transform=ax.transAxes,
fontsize=7 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
def DICOM_2D_image(ax, cmap, dicom_file, title, fontsize_adjustment=0, patien_age="__"):
ax.imshow(dicom_file.pixel_array, cmap)
ax.grid(False)
ax.set_title(title, loc="center", pad=10,
fontsize=12 + fontsize_adjustment)
ax.set(xticks=[], yticks=[], frame_on=False)
get_patient_information(ax, cmap, dicom_file,
fontsize_adjustment, patien_age)
get_plane_information(ax, cmap, dicom_file, fontsize_adjustment)
get_date_information(ax, cmap, dicom_file, fontsize_adjustment)
get_manufacturer_information(ax, cmap, dicom_file, fontsize_adjustment)
fig, axs = plt.subplots(2, 2, figsize=(
figure_size[0] / 3 * 2, figure_size[1] / 4 * 5))
DICOM_2D_image(
axs[0, 0],
"gray",
dicom_file,
f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'gray'",
)
DICOM_2D_image(
axs[0, 1],
"bone",
dicom_file,
f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'bone'",
)
DICOM_2D_image(
axs[1, 0],
"gray",
random_dicom_file,
f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'gray'",
)
DICOM_2D_image(
axs[1, 1],
"bone",
random_dicom_file,
f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'bone'",
)
fig.suptitle(
"Visualization of CT Medical Image Display with Supplementary Information",
fontsize="x-large",
x=0.5,
y=0,
)
plt.tight_layout()
plt.show()
# The dataset used in this practical project is the MRI DICOM dataset of the head of
# a 52-year-old normal male mathematics professor
# The subject of the experiment is the author, who suffers from small vertical strabismus
# (hypertropia), which is a misalignment of the eyes, which can be seen in this dataset
# This dataset was provided by Lionheart, William R.B. in 2015, available online
# (https://zenodo.org/record/16956 or http://doi.org/10.5281/zenodo.16956),
# license type is Attribution-ShareAlike 4.0 International (CC BY-SA 4.0)
# https://creativecommons.org/licenses/by-sa/4.0
path_to_head_mri = Path(
"../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
"ST000000/SE000002/"
)
# `Path.glob` selects the given relative pattern in the directory represented by the given path
# and yields all matching files (of any type)
# Since `Path.glob` returns a generator, it is converted to a list here for easy display
all_files = list(path_to_head_mri.glob("*"))
mri_data = []
for path in all_files:
# Read the individual DICOM files one by one in the list returned by the path generator
data = pydicom.read_file(path)
mri_data.append(data)
tabulation = Form_Generator()
tabulation.heading_printer("Reading all DICOM files that match the specified pattern")
statements = [
"""
path_to_head_mri = Path(
"../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
"ST000000/SE000001/"
)
all_files = list(path_to_head_mri.glob("*"))
mri_data = []
for path in all_files:
data = pydicom.read_file(path)
mri_data.append(data)
"""
]
tabulation.statement_generator(statements)
variables = ["str(path_to_head_mri)", "all_files"]
values = [str(path_to_head_mri), str(reprlib_rules.repr(all_files))]
tabulation.variable_generator(variables, values, 1)
expressions = [
"all_files[0].parts",
"all_files[0].parent",
"all_files[0].name",
"all_files[0].suffixes",
"all_files[0].stem",
"len(all_files)",
"mri_data[0]",
"len(mri_data)",
]
results = [
str(all_files[0].parts),
str(all_files[0].parent),
str(all_files[0].name),
str(all_files[0].suffixes),
str(all_files[0].stem),
str(len(all_files)),
str(reprlib_rules.repr(mri_data[0])),
str(len(mri_data)),
]
tabulation.expression_generator(expressions, results, 1)
Reading all DICOM files that match the specified pattern +---------------------------------------------------------+ | Statement | +---------------------------------------------------------+ | path_to_head_mri = Path( | | "../Datasets/An MRI Dicom Data Set of the Head of a | | Normal Male Human Aged 52/" | | "ST000000/SE000001/" | | ) | | all_files = list(path_to_head_mri.glob("*")) | | | | mri_data = [] | | | | for path in all_files: | | data = pydicom.read_file(path) | | mri_data.append(data) | +---------------------------------------------------------+ +-----------------------+-------------------------------------+ | Variable | Value | +-----------------------+-------------------------------------+ | str(path_to_head_mri) | ../Datasets/An MRI Dicom Data Set | | | of the Head of a Normal Male Human | | | Aged 52/ST000000/SE000002 | | all_files | [PosixPath('../Datasets/An MRI | | | Dicom Data Set of the Head of a | | | Normal Male Human Aged | | | 52/ST000000/SE000002/MR000000'), | | | PosixPath('../Datasets/An MRI | | | Dicom Data Set of the Head of a | | | Normal Male Human Aged | | | 52/ST000000/SE000002/MR000007'), | | | PosixPath('../Datasets/An MRI | | | Dicom Data Set of the Head of a | | | Normal Male Human Aged | | | 52/ST000000/SE000002/MR000031'), | | | PosixPath('../Datasets/An MRI | | | Dicom Data Set of the Head of a | | | Normal Male Human Aged | | | 52/ST000000/SE000002/MR000009'), | | | PosixPath('../Datasets/An MRI | | | Dicom Data Set of the Head of a | | | Normal Male Human Aged | | | 52/ST000000/SE000002/MR000008'), | | | PosixPath('../Datasets/An MRI | | | Dicom Data Set of the Head of a | | | Normal Male Human Aged | | | 52/ST000000/SE000002/MR000030'), | | | ...] | +-----------------------+-------------------------------------+ +-----------------------+-------------------------------------+ | Expression | Result | +-----------------------+-------------------------------------+ | all_files[0].parts | ('..', 'Datasets', 'An MRI Dicom | | | Data Set of the Head of a Normal | | | Male Human Aged 52', 'ST000000', | | | 'SE000002', 'MR000000') | | all_files[0].parent | ../Datasets/An MRI Dicom Data Set | | | of the Head of a Normal Male Human | | | Aged 52/ST000000/SE000002 | | all_files[0].name | MR000000 | | all_files[0].suffixes | [] | | all_files[0].stem | MR000000 | | len(all_files) | 32 | | mri_data[0] | Dataset.file_meta | | | ------------------------------- | | | (0002, 0000) File Meta Information | | | Group Length UL: 214 | | | (0002, 0001) Fil...Group Length | | | UL: 131084 | | | (7fe0, 0010) Pixel Data | | | OW: Array of 131072 | | | elements | | len(mri_data) | 32 | +-----------------------+-------------------------------------+
# DICOM files may not be sorted according to their actual image positions, this can
# be verified by checking Slice Location (0020,1041)
# Slice Location (0020,1041) identifies the relative position of the image plane in millimeters
unordered_slices = {
index: float(slice.SliceLocation) for index, slice in enumerate(mri_data)
}
# For better viewing and processing of 3D data stored with multiple 2D DICOM files,
# these files must be sorted; otherwise the complete scanned image will become cluttered
# and useless
mri_data_ordered = sorted(mri_data, key=lambda slice: slice.SliceLocation)
ordered_slices = {
index: float(slice.SliceLocation) for index, slice in enumerate(mri_data_ordered)
}
tabulation = Form_Generator()
tabulation.heading_printer("Sorting of read DICOM files")
statements = [
"""
unordered_slices = {
index: float(slice.SliceLocation) for index, slice in enumerate(mri_data)
}
mri_data_ordered = sorted(mri_data, key=lambda slice: slice.SliceLocation)
ordered_slices = {
index: float(slice.SliceLocation) for index, slice in enumerate(mri_data_ordered)
}
"""
]
tabulation.statement_generator(statements)
variables = ["unordered_slices", "ordered_slices"]
values = [
str(unordered_slices),
str(ordered_slices),
]
tabulation.variable_generator(variables, values, 1)
expressions = [
"mri_data_ordered[0]",
"len(mri_data_ordered)",
"len(unordered_slices)",
"len(ordered_slices)",
]
results = [
str(reprlib_rules.repr(mri_data_ordered[0])),
str(len(mri_data_ordered)),
str(len(unordered_slices)),
str(len(ordered_slices)),
]
tabulation.expression_generator(expressions, results, 1)
Sorting of read DICOM files +-----------------------------------------------------------+ | Statement | +-----------------------------------------------------------+ | unordered_slices = { | | index: float(slice.SliceLocation) for index, slice in | | enumerate(mri_data) | | } | | | | mri_data_ordered = sorted(mri_data, key=lambda slice: | | slice.SliceLocation) | | | | ordered_slices = { | | index: float(slice.SliceLocation) for index, slice in | | enumerate(mri_data_ordered) | | } | +-----------------------------------------------------------+ +------------------+------------------------------------------+ | Variable | Value | +------------------+------------------------------------------+ | unordered_slices | {0: 0.0, 1: 42.0000012462811, 2: | | | 185.999999047347, 3: 54.0000010703945, | | | 4: 48.000003062328, 5: | | | 180.000001085469, 6: 36.0000032341321, | | | 7: 5.99999800806651, 8: | | | 137.999999746811, 9: 144.000001609047, | | | 10: 72.0000008065647, 11: | | | 90.0000014947299, 12: 84.0000015826732, | | | 13: 78.0000016706165, 14: | | | 149.999999570924, 15: 132.000001784933, | | | 16: 23.9999996020383, 17: | | | 18.0000015939718, 18: 12.0000035859053, | | | 19: 30.0000014180852, 20: | | | 161.99999939912, 21: 108.000002136706, | | | 22: 120.00000196082, 23: | | | 96.0000013605981, 24: 174.000003031214, | | | 25: 66.000000894508, 26: | | | 60.0000009824512, 27: 102.000001319098, | | | 28: 168.000005069336, 29: | | | 126.000000015075, 30: 156.00000143316, | | | 31: 114.000002094951} | | ordered_slices | {0: 0.0, 1: 5.99999800806651, 2: | | | 12.0000035859053, 3: 18.0000015939718, | | | 4: 23.9999996020383, 5: | | | 30.0000014180852, 6: 36.0000032341321, | | | 7: 42.0000012462811, 8: | | | 48.000003062328, 9: 54.0000010703945, | | | 10: 60.0000009824512, 11: | | | 66.000000894508, 12: 72.0000008065647, | | | 13: 78.0000016706165, 14: | | | 84.0000015826732, 15: 90.0000014947299, | | | 16: 96.0000013605981, 17: | | | 102.000001319098, 18: 108.000002136706, | | | 19: 114.000002094951, 20: | | | 120.00000196082, 21: 126.000000015075, | | | 22: 132.000001784933, 23: | | | 137.999999746811, 24: 144.000001609047, | | | 25: 149.999999570924, 26: | | | 156.00000143316, 27: 161.99999939912, | | | 28: 168.000005069336, 29: | | | 174.000003031214, 30: 180.000001085469, | | | 31: 185.999999047347} | +------------------+------------------------------------------+ +-----------------------+-------------------------------------+ | Expression | Result | +-----------------------+-------------------------------------+ | mri_data_ordered[0] | Dataset.file_meta | | | ------------------------------- | | | (0002, 0000) File Meta Information | | | Group Length UL: 214 | | | (0002, 0001) Fil...Group Length | | | UL: 131084 | | | (7fe0, 0010) Pixel Data | | | OW: Array of 131072 | | | elements | | len(mri_data_ordered) | 32 | | len(unordered_slices) | 32 | | len(ordered_slices) | 32 | +-----------------------+-------------------------------------+
# Fill the 3D array in a slice-per-slice manner
full_volume = [slice.pixel_array for slice in mri_data_ordered]
full_volume = np.array(full_volume)
tabulation = Form_Generator()
tabulation.heading_printer(
"One-time extraction of the overall 3D pixel array for all slices"
)
statements = [
"""
full_volume = [slice.pixel_array for slice in mri_data_ordered]
full_volume = np.array(full_volume)
"""
]
tabulation.statement_generator(statements)
variables = ["full_volume"]
values = [str(reprlib_rules.repr(full_volume))]
tabulation.variable_generator(variables, values)
expressions = ["full_volume.shape", "full_volume.dtype"]
results = [str(full_volume.shape), str(full_volume.dtype)]
tabulation.expression_generator(expressions, results)
One-time extraction of the overall 3D pixel array for all slices +-----------------------------------------------+ | Statement | +-----------------------------------------------+ | full_volume = [slice.pixel_array for slice in | | mri_data_ordered] | | full_volume = np.array(full_volume) | +-----------------------------------------------+ +-------------+------------------------------------+ | Variable | Value | +-------------+------------------------------------+ | full_volume | array([[[0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | ..., | | | [0,... ..., | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0]]], | | | dtype=uint16) | +-------------+------------------------------------+ +-------------------+----------------+ | Expression | Result | +-------------------+----------------+ | full_volume.shape | (32, 256, 256) | | full_volume.dtype | uint16 | +-------------------+----------------+
def sort_checker(list_1):
list_2 = list_1[1:]
if all(a <= b for a, b in zip(list_1, list_2)):
return "ascending"
elif all(a >= b for a, b in zip(list_1, list_2)):
return "descending"
elif all(a == b for a, b in zip(list_1, list_2)):
return "identical"
else:
return "unordered"
def title_adder(func):
@wraps(func)
def wrapper(inputs_1, ax, dict_slices, title=None, **kwargs):
is_sorted = sort_checker(list(dict_slices.values()))
if title is None:
title = f"DICOM image grid with {is_sorted} slice locations"
func(inputs_1, ax, dict_slices, **kwargs)
ax.set_title(
title,
loc="center",
pad=10,
)
return wrapper
def grid_annotator(func):
@wraps(func)
def wrapper(inputs_1, ax, dict_slices, fontsize_adjustment=0, **kwargs):
x_positions_1, x_positions_2, y_positions, cmap = func(inputs_1, ax, **kwargs)
for (y, x_1), key in zip(
itertools.product(y_positions, x_positions_1), dict_slices.keys()
):
text = f"No.{(key + 1):>03}"
ax.text(
x_1,
y,
text,
transform=ax.transData,
fontsize=7 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
for (y, x_2), value in zip(
itertools.product(y_positions, x_positions_2), dict_slices.values()
):
text = f"{value:>5.1f} mm"
ax.text(
x_2,
y,
text,
transform=ax.transData,
fontsize=7 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
return wrapper
@title_adder
@grid_annotator
def grid_DICOM_2D_image(inputs_1, ax, n=None, row_size=9, pad_value=1):
if n is None:
n = inputs_1.shape[0]
nrows = math.ceil(n / row_size)
cmap = "bone"
# The only types supported by `torch.Tensor` are float64, float32, float16, complex64,
# complex128, int64, int32, int16, int8, uint8 and bool, not uint16
inputs_1 = torch.Tensor(inputs_1.astype(int)).unsqueeze(1)
x_positions_1 = [
math.ceil(
x * inputs_1[:row_size].shape[3]
+ 0.075 * inputs_1[:row_size].shape[3]
+ pad_value * x
)
for x in range(inputs_1[:row_size].shape[0])
]
x_positions_2 = [
math.ceil(
x * inputs_1[:row_size].shape[3]
+ 0.55 * inputs_1[:row_size].shape[3]
+ pad_value * x
)
for x in range(inputs_1[:row_size].shape[0])
]
y_positions = [
math.ceil(
y * inputs_1[:row_size].shape[2]
+ 0.125 * inputs_1[:row_size].shape[2]
+ pad_value * y
)
for y in range(nrows)
]
for i in range(nrows):
if len(inputs_1) > row_size:
row_images = inputs_1[:row_size]
inputs_1 = inputs_1[row_size:]
else:
row_images = inputs_1[:]
# Unlike the images in the MNIST dataset, the DICOM file images do not take values
# in the range (0, 1), so the parameter `normalize` needs to be set to True in order to
# move their values proportionally to the range (0, 1)
# By querying the source code of `vision/torchvision/utils.py`, it is clear that any
# single-channel image is converted to a three-channel image, and all are represented
# as tensors, as shown below:
# > if tensor.dim() == 2: # single image H x W
# > tensor = tensor.unsqueeze(0)
# > if tensor.dim() == 3: # single image
# > if tensor.size(0) == 1: # if single-channel, convert to 3-channel
# > tensor = torch.cat((tensor, tensor, tensor), 0)
# > tensor = tensor.unsqueeze(0)
# >
# > if tensor.dim() == 4 and tensor.size(1) == 1: # single-channel images
# > tensor = torch.cat((tensor, tensor, tensor), 1)
grid_row = make_grid(
row_images, nrow=row_size, normalize=True, pad_value=pad_value
)
scalar_row = grid_row.numpy().transpose(1, 2, 0)[:, :, 0]
if i == 0:
if nrows != 1:
scalar_image = scalar_row
else:
scalar_image = np.full(
(
scalar_row.shape[0],
math.ceil(scalar_row.shape[1] / inputs_1.shape[0] * row_size),
),
np.nan,
)
scalar_image[:, : scalar_row.shape[1]] = scalar_row
elif scalar_image.shape[1] == scalar_row.shape[1]:
scalar_image = np.concatenate((scalar_image, scalar_row), axis=0)
else:
image_to_concat = np.full(
(scalar_image.shape[0] + scalar_row.shape[0], scalar_image.shape[1]),
np.nan,
)
image_to_concat[: scalar_image.shape[0], :] = scalar_image
image_to_concat[scalar_image.shape[0] :, : scalar_row.shape[1]] = scalar_row
scalar_image = image_to_concat
ax.imshow(scalar_image, cmap=cmap)
ax.set(xticks=[], yticks=[], frame_on=False)
ax.grid(False)
return x_positions_1, x_positions_2, y_positions, cmap
unordered_full_volume = [slice.pixel_array for slice in mri_data]
unordered_full_volume = np.array(unordered_full_volume)
descending_slices = {
index: float(slice.SliceLocation)
for index, slice in enumerate(mri_data_ordered[::-1])
}
fig, axs = plt.subplots(3, 1, figsize=(figure_size[0], figure_size[1] / 7 * 9))
grid_DICOM_2D_image(unordered_full_volume, axs[0], unordered_slices)
grid_DICOM_2D_image(full_volume, axs[1], ordered_slices)
grid_DICOM_2D_image(full_volume[::-1], axs[2], descending_slices)
fig.suptitle(
"Visual Comparison of DICOM Image Grids with Unordered and Ordered Slice Locations",
fontsize="x-large",
x=0.5,
y=0,
)
plt.tight_layout()
plt.show()
fig, axs = plt.subplots(9, 3, figsize=(
figure_size[0] / 3 * 2, figure_size[1] / 2 * 7))
for slice_counter, (i, j) in enumerate(itertools.product(range(9), range(3))):
DICOM_2D_image(
axs[i][j],
"bone",
mri_data_ordered[slice_counter],
f"DICOM image with {mri_data_ordered[slice_counter].SliceLocation:^.1f} mm "
"slice location",
fontsize_adjustment=-3,
patien_age="52",
)
fig.suptitle(
"Visualization of MRI Medical Images with Information Displayed in Order of Slice Location",
fontsize="x-large",
x=0.5,
y=0,
)
plt.tight_layout()
plt.show()
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) Cell In[20], line 5 1 fig, axs = plt.subplots(9, 3, figsize=( 2 figure_size[0] / 3 * 2, figure_size[1] / 2 * 7)) 4 for slice_counter, (i, j) in enumerate(itertools.product(range(9), range(3))): ----> 5 DICOM_2D_image( 6 axs[i][j], 7 "bone", 8 mri_data_ordered[slice_counter], 9 f"DICOM image with {mri_data_ordered[slice_counter].SliceLocation:^.1f} mm " 10 "slice location", 11 fontsize_adjustment=-3, 12 patien_age="52", 13 ) 15 fig.suptitle( 16 "Visualization of MRI Medical Images with Information Displayed in Order of Slice Location", 17 fontsize="x-large", 18 x=0.5, 19 y=0, 20 ) 22 plt.tight_layout() Cell In[15], line 264, in DICOM_2D_image(ax, cmap, dicom_file, title, fontsize_adjustment, patien_age) 261 ax.set(xticks=[], yticks=[], frame_on=False) 262 get_patient_information(ax, cmap, dicom_file, 263 fontsize_adjustment, patien_age) --> 264 get_plane_information(ax, cmap, dicom_file, fontsize_adjustment) 265 get_date_information(ax, cmap, dicom_file, fontsize_adjustment) 266 get_manufacturer_information(ax, cmap, dicom_file, fontsize_adjustment) Cell In[15], line 30, in plane_information_annotator.<locals>.wrapper(ax, cmap, dicom_file, fontsize_adjustment) 28 @wraps(func) 29 def wrapper(ax, cmap, dicom_file, fontsize_adjustment): ---> 30 plane = func(ax, cmap, dicom_file, fontsize_adjustment) 31 text = f"Image: {dicom_file.Columns} x {dicom_file.Rows}" 32 text += f"\nImage Plane: {plane}" Cell In[15], line 52, in orientation_annotator.<locals>.wrapper(ax, cmap, dicom_file, fontsize_adjustment) 50 @wraps(func) 51 def wrapper(ax, cmap, dicom_file, fontsize_adjustment): ---> 52 plane, top, right, bottom, left = func(dicom_file) 53 ax.text( 54 0.4875, 55 0.95, (...) 61 c=annotation_color(cmap), 62 ) 63 ax.text( 64 0.4875, 65 0.025, (...) 71 c=annotation_color(cmap), 72 ) Cell In[15], line 143, in image_orientation.<locals>.wrapper(*args, **kwargs) 141 if patient_position.endswith("S"): 142 top = 0 --> 143 right = term_list.index(patient_position[:2]) * 2 + 1 144 left = right + 2 145 else: ValueError: 'HF' is not in list
# The path object must be converted to a string in order for the SimpleITK package to process
# the path information
path_string = str(path_to_head_mri)
# For some image formats such as DICOM, the images also contain associated metadata,
# which is not loaded by the reader by default for time-saving reasons
# `ImageSeriesReader` class provides the ability to read a series of image files into a
# SimpleITK image, and once a series of images has been read, the metadata can be accessed
# directly from the reader
# `GetGDCMSeriesIDs` provides the ability to get all seriesIDs from a DICOM dataset
series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(path_string)
# `GetGDCMSeriesFileNames` generates a series of file names from the directory containing
# the DICOM dataset and series ID
# Note that the resulting series of file names (which are actually file paths) are
# automatically sorted
series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(
path_string, series_ids[0]
)
series_reader = sitk.ImageSeriesReader()
series_reader.SetFileNames(series_file_names)
# Finally, the reader is executed by calling `Execute` to get the required data
image_data = series_reader.Execute()
tabulation = Form_Generator()
tabulation.heading_printer(
"Automatic recognition and recall of series file names for DICOM images"
)
statements = [
"""
path_string = str(path_to_head_mri)
series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(path_string)
series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(
path_string, series_ids[0]
)
series_reader = sitk.ImageSeriesReader()
series_reader.SetFileNames(series_file_names)
image_data = series_reader.Execute()
"""
]
tabulation.statement_generator(statements)
variables = ["series_ids", "series_file_names", "series_reader", "image_data"]
values = [
str(reprlib_rules.repr(series_ids)),
str(reprlib_rules.repr(series_file_names)),
str(reprlib_rules.repr(series_reader)),
str(reprlib_rules.repr(image_data)),
]
tabulation.variable_generator(variables, values, 1)
expressions = [
"type(series_ids)",
"len(series_ids)",
"type(series_file_names)",
"len(series_file_names)",
"len(image_data)",
"image_data.GetSize()",
]
results = [
str(type(series_ids)),
str(len(series_ids)),
str(type(series_file_names)),
str(len(series_file_names)),
str(len(image_data)),
str(image_data.GetSize()),
]
tabulation.expression_generator(expressions, results)
Automatic recognition and recall of series file names for DICOM images +----------------------------------------------------------+ | Statement | +----------------------------------------------------------+ | path_string = str(path_to_head_mri) | | series_ids = | | sitk.ImageSeriesReader.GetGDCMSeriesIDs(path_string) | | series_file_names = | | sitk.ImageSeriesReader.GetGDCMSeriesFileNames( | | path_string, series_ids[0] | | ) | | | | series_reader = sitk.ImageSeriesReader() | | series_reader.SetFileNames(series_file_names) | | | | image_data = series_reader.Execute() | +----------------------------------------------------------+ +-------------------+-----------------------------------------+ | Variable | Value | +-------------------+-----------------------------------------+ | series_ids | ('1.3.46.67058...1413262801702',) | | series_file_names | ('../Datasets/...0001/MR000000', | | | '../Datasets/...0001/MR000001', | | | '../Datasets/...0001/MR000002', | | | '../Datasets/...0001/MR000003', | | | '../Datasets/...0001/MR000004', | | | '../Datasets/...0001/MR000005', ...) | | series_reader | <SimpleITK.SimpleITK.ImageSeriesReader; | | | proxy of <Swig Object of type | | | 'itk::simple::ImageSeriesReader *' at | | | 0x15a2ff3f0> > | | image_data | <SimpleITK.SimpleITK.Image; proxy of | | | <Swig Object of type 'std::vector< | | | itk::simple::Image >::value_type *' at | | | 0x15a2fced0> > | +-------------------+-----------------------------------------+ +-------------------------+-----------------+ | Expression | Result | +-------------------------+-----------------+ | type(series_ids) | ⟨class 'tuple'⟩ | | len(series_ids) | 1 | | type(series_file_names) | ⟨class 'tuple'⟩ | | len(series_file_names) | 27 | | len(image_data) | 1769472 | | image_data.GetSize() | (256, 256, 27) | +-------------------------+-----------------+
# As you can see above, the original SimpleITK image has the shape (256, 256, 27) instead of
# the common shape (27, 256, 256), this is due to the different order of the image sizes,
# so it is necessary to perform the last step to convert the SimpleITK image to a NumPy array
# `GetArrayFromImage` returns a copy of the image data, which can be freely modified
# as it has no effect on the original SimpleITK image
head_mri = sitk.GetArrayFromImage(image_data)
tabulation = Form_Generator()
tabulation.heading_printer("Getting the NumPy array from the SimpleITK image")
statements = ["head_mri = sitk.GetArrayFromImage(image_data)"]
tabulation.statement_generator(statements)
variables = ["head_mri"]
values = [str(reprlib_rules.repr(head_mri))]
tabulation.variable_generator(variables, values)
expressions = [
"type(head_mri)",
"len(head_mri)",
"head_mri.shape",
]
results = [
str(type(head_mri)),
str(len(head_mri)),
str(head_mri.shape),
]
tabulation.expression_generator(expressions, results)
Getting the NumPy array from the SimpleITK image +-----------------------------------------------+ | Statement | +-----------------------------------------------+ | head_mri = sitk.GetArrayFromImage(image_data) | +-----------------------------------------------+ +----------+--------------------------------------------------+ | Variable | Value | +----------+--------------------------------------------------+ | head_mri | array([[[0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | ..., | | | [0,... ..., | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0]]], dtype=uint16) | +----------+--------------------------------------------------+ +----------------+-------------------------+ | Expression | Result | +----------------+-------------------------+ | type(head_mri) | ⟨class 'numpy.ndarray'⟩ | | len(head_mri) | 27 | | head_mri.shape | (27, 256, 256) | +----------------+-------------------------+
head_mri_directory_path = Path(
"../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
"ST000000/"
)
all_paths = sorted(list(head_mri_directory_path.glob("SE*")))
path = all_paths[-1]
new_series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(str(path))
new_series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(
str(path), new_series_ids[0]
)
new_series_reader = sitk.ImageSeriesReader()
new_series_reader.SetFileNames(new_series_file_names)
new_image_data = new_series_reader.Execute()
new_head_mri = sitk.GetArrayFromImage(new_image_data)
new_all_files = list(path.glob("*"))
new_mri_data = sorted(
[pydicom.read_file(p) for p in new_all_files],
key=lambda slice: slice.SliceLocation,
)
new_ordered_slices = {
index: float(slice.SliceLocation) for index, slice in enumerate(new_mri_data)
}
fig = plt.figure(figsize=(figure_size[0], figure_size[1]), constrained_layout=True)
gs = gridspec.GridSpec(
nrows=2, ncols=1, figure=fig, wspace=0.08, hspace=0.02, height_ratios=[3, 4]
)
ax_1 = fig.add_subplot(gs[0, 0])
grid_DICOM_2D_image(head_mri, ax_1, ordered_slices, title="xxx")
ax_2 = fig.add_subplot(gs[1, 0])
grid_DICOM_2D_image(new_head_mri, ax_2, new_ordered_slices, title="xxx")
fig.suptitle(
"xxx",
fontsize="x-large",
x=0.5,
y=0,
)
plt.show()